%% it solves the LCP problem

BB                = sparse((r+delta)*eye(N)- diag(drift_I,0))-AA;
bb                = uu;
S                 = v_tresh;
[w,s_new,retcode] = LCPSolve(BB,-bb+BB*S);
LCP_error     = s_new.*(BB*s_new -bb+BB*S);
LCP_error_max = max(abs(LCP_error));
if LCP_error_max > 1/100000
    disp('LCP not solved')
    return
end
v_lcp      = s_new+S;